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A bstract 

The  purpose  of  this  study  is  to  locate  critical  inclinations  in  long  term  high  eccen¬ 
tricity  orbits  about  Mars  using  numerical  methods.  A  critical  inclination  is  defined  as  the 
inclination  at  orbit  insertion  which  produces  a  local  maximum  in  the  amplitude  of  the  vari¬ 
ation  of  eccentricity  or  inclination.  The  perturbation  model  consists  of  the  first  non-zero 
zonal  harmonic  ( J2)  of  the  Mars  gravity  potential  and  the  Sun  as  a  point  mass  third  body. 
The  search  range  consists  of  inclinations  from  0.25  to  90.0  degrees,  eccentricities  from  0.40 
to  0.90,  periapse  radii  from  4000  km  to  7000  km,  and  orbit  lifetimes  of  10  Earth  years. 

The  numerical  search  comprises  the  following  procedure;  (1)  A  time  history  of  eccen¬ 
tricity  and  inclination  is  produced  for  each  combination  of  orbit  insertion  initial  conditions 
by  numerically  propagating  Lagrange’s  Planetary  Equations.  (2)  Each  time  history  is  fit, 
in  the  least  squares  sense,  to  a  linear  function.  The  standard  deviation  of  the  residuals 
for  each  fit  is  employed  as  the  search  parameter.  (3)  A  three-dimensional  surface  plot 
of  the  standard  deviation  in  eccentricity  and  the  standard  deviation  in  inclination  versus 
eccentricity  and  inclination  is  produced  for  each  value  of  periapse  radius  considered.  The 
local  maximums  in  these  surfaces  identifies  the  locations  of  the  critical  inclinations.  (4) 
The  three-dimensional  surfaces  are  then  reduced  to  two  dimensions  by  plotting  inclination 
versus  eccentricity  for  the  local  maximums  in  stajidard  deviation. 

Six  critical  inclination  curves  for  eccentricity  are  identified,  three  of  which  are  curve 
fit  and  found  to  be  linear  in  periapse  radius  and  quadratic  in  eccentricity.  The  surface 
plots  for  inclination  indicate  the  presence  of  large  variations  but  the  surface  topography 
does  not  allow  for  the  identification  of  distinct  local  maximum  curves. 


NUMERICAL  DETERMINATION  OF  THE  LOCATION  OF  CRITICAL 
INCLINATIONS  FOR  LONG  TERM  HIGH  ECCENTRICITY 
ORBITS  ABOUT  MARS 


/.  Introduction 

The  task  of  placing  satellites  in  long-term  orbit  about  Mars  will  require  extensive 
planning,  but  one  of  the  first  mission  parameters  to  consider  will  be  the  orbit  itself.  Orbit 
design  will  be  a  function  not  only  of  the  specific  mission  requirements,  but  also  of  mission 
costs,  and  in  particular,  the  costs  related  to  orbit  insertion.  Inserting  a  satellite  into  orbit 
about  Mars  requires  a  large  velocity  change  between  the  Earth-to-Mars  transfer  trajectory 
and  a  closed  orbit  about  Mars.  The  magnitude  of  the  velocity  change,  and  therefore  the 
related  costs,  is  inversely  proportional  to  the  energy  of  the  Mars-centered  orbit. 

Consider  a  satellite  mission  which  includes  the  need  to  make  close  observations  of 
the  Martian  surface.  This  requirement  may  be  met  by  simply  inserting  the  satellite  into 
a  near  circular  orbit  just  above  the  atmosphere,  but  such  an  orbit  possesses  at  least  two 
disadvantages:  1)  atmospheric  drag,  and  2)  the  fuel  required  to  insert  the  satellite  into 
such  a  low-energy  orbit. 

An  approach  to  the  solution  of  these  problems  is  to  utilize  a  high  eccentricity  orbit, 
which  allows  for  close  approach  (at  periapse)  as  well  as  reduced  insertion  fuel  requirements 
(compared  to  the  circular  or  near  circular  orbit).  High  eccentricity  orbits,  however,  are 
not  without  difficulties.  The  combined  effects  of  the  gravity  fields  of  Mars  and  of  the  Sun 
will  produce  resonance,  causing  large  variations  in  the  orbital  elements.  Such  variations 
may  result  in  violated  mission  constraints,  and  in  severe  cases,  reentry  and  impact. 

To  illustrate  this  phenomenon,  Figures  1  and  2  display  the  behavior  of  the  eccentricity 
and  inclination  of  a  Martian  orbit  over  a  period  of  twenty  Earth  years  (under  the  conditions 
for  which  large  variations  are  not  present).  Figures  3  and  4  present  the  behavior  of  these 
parameters  under  the  conditions  for  which  large  variations  do  exist. 
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Figure  1.  Example  of  Small  Amplitude  Variation  in  Eccentricity  for  a  Long-term  Martian 
Orbit 
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Figure  2.  Example  of  Small  Amplitude  Variation  in  Inclination  for  a  Long-term  Martian 
Orbit 
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Example  of  Large  Amplitude  Variation  in  Eccentricity  for  a  Long-term  Martian 
Orbit 
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Example  of  Large  Amplitude  Variation  in  Inclination  for  a  Long-term  Martian 
Orbit 


The  conditions  for  which  such  large  variations  occur  are  primarily  a  function  of  the 
values  of  the  eccentricity  and  inclination  at  orbit  insertion.  These  two  elements  determine 
the  exposure  of  the  orbit  to  the  effects  of  the  planet  oblateness  and  to  the  third-body 
effects  of  the  Sun.  The  dependence  of  these  variations  on  ectentricity  will  be  shown  to 
be  generally  a  direct  one,  that  is,  the  larger  the  eccentricity,  the  greater  the  variations. 
However,  the  effects  of  the  orbit  inclination  are  not  as  clear,  and  are  therefore  of  greater 
concern. 

This  thesis  develops  a  numerical  method  for  determining  the  values  of  inclination  and 
eccentricity  for  which  large  variations  in  these  two  elements  occur.  The  gravity  models 
will  be  simplified,  considering  only  the  first  non-zero  zonal  harmonic  (J2  )  of  the  Mars 
gravity  potential,  and  treating  the  Sun  as  a  point  mass.  The  reasons  for  this  are  twofold: 
1)  to  meet  with  the  constraint  of  limited  computer  time,  and  2),  to  provide  a  basis  for 
comparison  with  an  analytical  approach  in  locating  the  critical  inclinations.  This  model 
therefore  includes  only  the  dominant  sources  of  perturbation  (for  Mars,  J2  is  larger  than 
any  other  zonal  harmonic  by  two  orders  of  magnitude,  and  is  twice  the  magnitude  of 
for  Earth). 

This  investigation  begins  by  presenting  the  theory  upon  which  the  numerical  work 
is  based,  including  the  derivation  of  the  disturbing  functions  used  in  Lagrange’s  Planetary 
Equations  (LPEs)  for  propagating  the  time  history  of  the  orbital  elements,  as  well  as 
providing  a  brief  summary  of  an  analytical  approach  to  locating  the  critical  inclinations. 
The  numerical  work  is  presented  by  describing  the  computer  program  which  was  utilized, 
followed  by  a  description  of  the  approach  in  locating  the  critical  inclinations.  The  results 
are  given  primarily  in  graphical  form. 


II.  Analytical  Development 


Equation  of  Motion 

Figure  5  shows  the  reference  frame  used  throughout  this  investigation.  The  Sun, 
relative  to  this  reference  frame,  is  considered  to  “orbit”  Mars. 


Figure  5.  Mars-centered  Reference  Frame 
Symbols  used  in  this  figure  are  defined  as  follows: 


r  =  satellite  radius  vector 
r,  =  Sun  radius  vector 

9  i  =  inclination  of  satellite  orbit 

f2  =  argument  of  ascending  node  for  satellite  orbit 
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u;  =  argument  of  periapse  for  satellite  orbit 
I  =  inclination  of  Sun  orbit 
A  =  right  ascension  of  the  Sun 
L  =  Mars-centered  latitude  of  the  Sun 
6  =  satellite  central  angle 
6g  =  Sun  central  angle 

The  equations  of  motion  for  the  satellite  in  this  reference  frame  may  be  written  as 
the  gradient  of  the  gravity  potentials  of  Mars  and  of  the  Sun: 

f  =  nRM+Rs)  (1) 

where 

V  =  gradient  operator 
Rm  =  Mars  gravity  potential 
Rs  =  Sun  gravity  potential 

The  orbital  elements  (other  than  the  mean  anomaly)  remain  constant  when  the 
central  body  is  modeled  as  a  point  mass  and  no  other  forces  are  included.  The  effects 
of  the  S\m  and  of  the  oblateness  of  Mars  introduce  perturbations  which  cause  the  orbital 
elements  to  vary  with  time.  The  perturbing  contribution  to  the  gravity  potentials  are 
modeled  as  a  disturbing  function  R  ,  which,  when  introduced  into  the  LPEs  [4:476-483], 
describes  the  ‘motion’  of  the  orbital  elements. 


The  LPEs  may  be  written  as  follows: 


dn 

1  dR 

dt 

nab  sin  i  di 

(2) 

di 

1  dR  _  cost  dR 

dt 

= 

nab  sin :  dil  nab  sin  t  du 

(3) 

du 

cost  dR  _  6  dR 

dt 

nabsini  di  na^e  de 

(4) 

da 

2  dR 

dt 

na  dM 

(5) 

de 

b^  dR  b  dR 

dt 

na*e  dM  na^e  du 

(6) 

dM 

2  dR  6*  dR 

dt 

* 

na  da  na*e  de 

(7) 

where 

e  =  eccentricity  of  satellite  orbit 

a  =  semi-major  axis  of  satellite  orbit 

b  =  ay/1  -  =  semi-minor  ajds  of  satellite  orbit 

Mm  =  Mars  gravity  constant 

M  =  satellite  mean  anomaly 

n  =  =  satellite  mean  motion 

The  disturbing  function  Ji  will  now  be  derived. 


7 


The  Mars  Disturbing  Function 


The  gravity  potential  of  Mars  may  be  written  as  a  combination  of  the  two-body, 
undisturbed  potential  and  the  disturbing  function: 

Rm  =  ~F  Rm  (8) 

T 

where 


T  =  magnitude  of  r 

Rm  =  Mars  disturbing  function 


The  Mars  disturbing  function  (3:421),  including  only  the  second  harmonic  coefficient  J2 
(i.e.,  the  first  non-zero  zonal  harmonic),  may  be  written  as 


Rm  = 


f^mJ2Rl 


(1  —  3sin^  6) 


where  Rg  is  the  Mars  equatorial  radius,  and  S  is  the  satellite  declination  (see  Figure  5), 


The  parameter  S  must  now  be  rewritten  in  terms  of  the  orbital  elements  /  and  w, 
where  /  is  the  satellite  true  anomaly.  To  carry  out  this  transformation,  consider  Figures  6 
and  7  on  the  following  page. 


Figure  6.  Looking  Down  the  Mars  Spin  Axis 


Figure  7.  Relationship  between  6  aoid  f  -j-u  :  A  Right  Spherical  Triangle 


In  Figure  6,  u;  is  measured  counterclockwise  from  the  ascending  node,  /  is  measured 
counterclockwise  from  periapse,  and  a  is  defined  such  that  a  +  /  +  w  =  tt.  The  central 
angles  a  and  6  and  the  vertex  angle  t  define  a  right  spherical  triangle,  as  shown  in  Figure  7. 
Making  use  of  spherical  and  plane  trigonometric  identities  gives 

sin  S  =  sin  t  sin  a 

=  sintsin(7r  —  (/ +  w)) 

=  sin  i  sin(/  +  w) 

Squaring  yields 

sin^  5  =  sin^  t  sin^(/ +  a>) 

=  sin^  i  |l  —  cos^(/  +  u;)j 
=  sin^t  [l  -  Q  +  icos2(/  +  u;)^j 

=  i  sin^  t  [1  -  cos  2(/  +  w)] 

Substituting  for  sin^  6  in  equation  (9)  results  in 

~  ^  ^  I  *  cos2(/  +  u>)^  (10) 


Averaging  the  Mars  Disturbing  Function 

To  determine  the  long-term  behavior  of  the  orbital  elements,  the  short-term  effects 
of  the  mean  (or  true)  anomaly  may  be  averaged  out  by  integrating  the  disturbing  function 
over  2-k  with  respect  to  the  mean  cinomaly  M.  The  averaged  Mars  disturbing  function  Ryn 
is  acquired  by 

r  RmdM  (11) 

Jo 

This  integration  is  fascillitated  by  employing  Hansen’s  Coefficients  [2].  A  brief  description 
of  Hansen’s  Coefficients  along  with  the  derivation  of  the  Hansen’s  Coefficients  used  below 
are  given  in  Appendix  A.  Substituting  equation  (10)  into  equation  (11)  yields 


2jr  Jo  2r3  V  2  J 


R. 


dM 


+ 


27r  Jo  i 


2r3 

2^  ijlrnJiRl  ,•  2(/  -b  W)  dM 

4rJ 


Calling  the  first  integral  A  and  the  second  B: 

Hn^hRl 


A  = 


2a3 


(12) 


(13) 


From  appendix  A: 


Therefore 


-3/2 


A  = 


2a3(l-e2)^/* 
Working  with  the  second  integral: 

Zh^J2RI 


(14) 


(15) 


B  = 


4a3 


“*2(/  +  a,)dM|  (16) 


Using  the  appropriate  trigonometric  relation,  the  expression  in  curly  brackets  becomes 
(0  (cos2/cos2u;  -  sin2/sin2u»)  dM 

‘  r  iz)  ^  r  (0 '  ““  " 


=  cos; 
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From  appendix  A 


and 


1  /■25r  /  ^  \  -3 

—  /  (  —  I  sin2/dAf  =  0 

27r  Jo  \a  J 


Hence,  the  averaged  Mars  disturbing  function  becomes 

-TT-  timhRl  ( ^  3  •  2  A 

2a3(l V  2  } 


(17) 

(18) 

(19) 
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The  Solar  Disturbing  Function 


Since  the  third-body  effects  of  the  Sun  are  to  be  included  by  considering  the  Sun  to 
be  a  point  mass,  the  solar  disturbing  function  may  be  derived  by  invoking  Newton’s  Law 
of  Universal  Gravitation.  The  final  expression  for  the  solar  disturbing  function  is  to  be 
with  repect  to  the  reference  frame  shown  in  Figure  5.  To  accomplish  this,  consider  the 
arbitrary  inertial  reference  frame  shown  in  Figure  8  below. 


z 


Figure  8.  Inertial  Reference  Frame  for  Deriving  Solar  Disturbing  Function 
Define  the  following: 

fi  =  radius  vector  to  satellite 
fi  =  radius  vector  to  Mars 
fa  =:  radius  vector  to  Sun 
mi  =  mass  of  satellite 


The  force  on  the  satellite  due  to  the  combined  effects  of  Mars  and  the  Sun  is 


Fi  =  —Gm\ 


—T21  +  -3-r3i 


(20) 


'21  '31 

where  G  is  the  universal  gravitational  constant.  The  force  on  Mars  due  to  the  combined 
effects  of  the  satellite  and  the  Sun  is 


F2  ~  — Gttij 


Wii_  ,  m3_ 

+  rr’‘32 
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'32 


The  accelerations  of  the  satellite  and  of  Mars  are 


mi 

f2  =  —  =  -G 
m2 


m2  m3 

^T21  +  ^3^731 


'21 


'31 


m2  m3 

.  T‘^12  +  “3"”32 

Lri2  r^2 


(21) 


(22) 


(23) 


Since  r2i  =  ri  -  r2,  the  acceleration  of  the  satellite  relative  to  Mars  is  7-21  =  fi  -  7*2,  or 


^21  =  -G 


m2^  m3_ 

Zr^2i  +  -r^3i 


'21 


'31 


+  G 


fmi  m3 

^ri2  +  ^r32 

1^12  C32 


=  -G 


(m2  +  mi) 


F21  +  Gms 


'21 


T32  rsi 


'32 


(24) 


Since  m2  >  mi,  mi  may  be  ignored.  Making  the  substitutions  fim  =  Gm2  and  /i,  =  Gmz 
results  in 

r2i  =  -^r2i-/i,  N--p-  (25) 

7'21  17^31  7^32  J 

In  terms  of  the  notation  used  in  Figure  5 


f'21 


T 


T23  =  -1^32  = 


rsi  =  r  -  r. 

Making  these  substitutions,  equation  (25)  becomes 


^  =  lT=-^  +  — 


T-fs 


'  r? 


=  ^{Rm  +  Rs) 


Therefore 


_  r  7,-7  7s 
rl 


(26) 


(27) 
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The  solar  gravity  potential  is  then  (see  (11)) 

r,  .  r  1 

Rs  —  I  3 

\.\rs-r\  rf 


(28) 


The  solar  disturbing  function  iZ,  may  now  be  extracted  from  the  solar  gravity 
potential  Rs-  Let  p  =  r,  ~r.  Therefore 

=  |r,  _  rp  =  (r,  -  r)  •  (r,  -  r) 

=  r  j  •  Tj  —  2r  •  r j  +  r  •  r 
=  ^  2rrs  cos  B  + 


where 


Factoring  out 


cosB  = 


T  •  r. 


rr. 


(29) 


2  2  /  ,  2r  ' 

=  r,  (  1+  - cos  5 


P  =  T» 


T1/2 


1  +  -r - COS  B 


=  I?,  ~  ?| 


(30) 


Substituting  Eqs  (29)  and  (30)  into  Eq  (28)  yields 


T. 


\  2r 

1  +  -=■ - cos  B 


-1/2 


rcos5 


(31) 


The  first  term  within  the  brackets  of  Eq  (31)  may  be  expanded  in  a  binomial  series: 


1  /r*  2r  \  1 3  /r*  2r  \ 

=  1-:;  -0-— cos  5  cos5  + 


2  V  r2  r. 


2  4  \  r2  r. 


Since  r,  >  r,  terms  of  order  (r/r,)^  and  above  may  be  neglected,  yielding 


1  +  -7 - cos  B 


-1/2 


1  /r*  2r  \  3/4r2  \ 


(32) 
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Substituting  Eq  (32)  into  Eq  (31)  results  in 


Rs  =  — 

r. 


Ms 


r‘  3r2 


l  +  ^(3co»“B-l) 


Analogous  to  Eq  (8) 

Therefore,  the  solar  disturbing  function  is 


Rs  =  —  +  Rs 

T, 


-  2rJ 


or 


|3  cos^  5  —  ij 

ra'- 


(33) 


(34) 


(35) 


(36) 


The  parenthetical  expression  in  Eq  (36)  must  now  be  written  in  terms  of  the  elements 
of  the  satellite  orbit  and  of  the  apparent  Sun  orbit.  In  order  to  simplify  notation,  the 
following  abbreviations  will  be  used  henceforth  for  all  trigonometric  functions: 

sin  0  ~  Se  cos  0  =  Cg 

Using  the  notation  of  Figure  5,  define  the  following  unit  vectors: 


€r  — '  ' 


CfiCg  -  SfiCiSg 
SnCg  +  CnCiSg 
SiSg 

C\Ci 
ClS\ 

Sl 


(37) 


(38) 


Therefore 


r  •  r. 


rr. 


—  Cr  •  Cm 


=  C^CgC^Ci,  —  SnCiSgC^Ci  +  SfaCgCiS/^^  +  CnCiSgCiS^  +  SiSgSi 
=  CgCi  (CqCa  +  SuSa)  +  CiSgCi  {C^Sa  -  SqCa)  +  SiSgSi, 
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=  QC'i/C'n-A  +  CiSeCiSci-A  +  SiSeSi 


The  solar  disturbing  function  now  becomes 


2 

^  {3  [CeCiCci.^  +  CiS»CLSn.A  +  S.SeSi,]^  -  l} 


Expanding  the  term  in  brackets  yields 


where 


2 

=  [KiC]  +  2K2CeSe  +  KsSj]  -  l} 


=  (ClCu^a) 

=  SiCiSiCa^ti.  -  ChCiSii^ACu-A 
=  {CiCiSu-A  -  SLSif 


j  +  5<7w 


i-ic« 


2SgCg  — 


9  =  f  +  u; 


so  the  solar  disturbing  function  is  finally 


■R*  =  +  3>l252(/+u,)  +  2^3  “ 


where 


Ai  =  Ki-Kz 
A2  =  K2 


A3  =  K\Ar  Kz 
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Averaging  the  Solar  Disturbing  Function 

As  was  done  in  the  case  of  the  Mars  disturbing  function,  the  short-term  effects  of  the 
satellite  true  anomaly  may  be  averaged  out.  The  averaged  solar  disturbing  function  is 
obtained  by 

r  RsdM  (43) 

Jo 

To  accomplish  the  averaging  by  employing  Hansen’s  Coefficients  (appendix  A),  the  true 
anomaly  /  must  be  separated  from  the  argument  of  periapse  u.  Using  the  trigonometric 
relations 


C'2(/+w)  =  —  S2jS2w 

•^2(/+w)  =  ■‘r  C2j  Siu 


Eq  (42)  becomes 

R,  =  I -AiC2/C2a;  -  -v4i52/52u.  +  3A252/C'2u/ +  342^2/521^  + -/is  -  l| 

=  I  3A252u.^  C2/  A  ^3A2C’2u;  -  -Ai52u,^  S^j  +  -  A3  -  l| 

Multiplying  and  dividing  by  the  square  of  the  semi-major  axis  o,  and  then  distributing  r* 
results  in 


E,  = 


/x,a‘ 

2r3 


-A1C2U1  +  3A2S2, 


(= 


-t-  I  3A2C2W  -  -Ai52w 
Averaging  the  terms  which  are  functions  of  the  mean  anomaly  yields 

f2-r  /^  \  2  _  _  t;^2 


The  averaged  solar  disturbing  function  is  then 


|-^e^AiC2u;  +  15e3A252u>  +  (2  -|-  3e^^  ^2“^^  ~ 


4r3 


(44) 


(45) 

(46) 

(47) 


(48) 
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Substituting  for  Aj,  A2,  and  A3  yields 


17  = 


^  { ye"  [Clch_j,  -  (CiCLSn-A  -  5x,5if ]  Cj,. 

+  15e^  [CiS^SjCo-A  -  C'i<^i<^n-A'S'n-A]  •^20; 

+  (2  +  3e2)  [I  [ciC^A  +  -  SiSif]  -  l]  } 


(49) 


Eqs  (19)  and  (49)  combine  to  form  the  complete  disturbing  function  being  considered 
in  this  investigation: 

iZ  =  177  +  17  (50) 
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Analytical  Determination  of  the  Location  of  Critical  Inclinations 

The  last  section  detailed  the  derivation  of  the  disturbing  function  which  will  be  in¬ 
corporated  into  the  LPEs  for  numerical  integration.  The  resulting  disturbing  function  was 
obtained  by  averaging  on  2x  over  the  satellite  mean  anomaly  for  both  the  Mars  disturbing 
function  and  the  Solar  disturbing  function.  What  remains,  however,  are  the  time  depen¬ 
dent  terms  of  the  Sun’s  motion,  namely  the  solar  central  angle  9^  and  the  solar  latitude 
L  (see  Figure  5).  Averaging  with  respect  to  these  terms  fascillitates  the  development  of 
an  analytical  approach  for  determining  the  location  of  the  critical  inclinations  (see  (5:182- 
192)).  The  purpose  of  this  section  is  to  examine  the  conditions  under  which  such  averaging 
is  valid,  and  to  summarize  the  results  of  this  approach. 

The  averaging  of  any  cyclic  parameter  is  valid  when  its  variation  is  faster  than  the 
variation  of  the  other  variables  involved.  Averaging  with  respect  to  the  mean  anomaly  is 
valid  since,  in  this  investigation,  the  orbit  period  is  shorter  than  the  period  of  the  variation 
of  the  argument  of  periapse  w  or  of  the  argument  of  the  ascending  node  ft.  To  average 
out  the  effects  of  0,  and  £,  the  mean  motion  of  the  Sun  n,  must  also  be  greater  than  the 
rate  at  which  w  and  fl  vary.  The  conditions  for  which  this  relationship  holds  will  now  be 
demonstrated. 


To  obtain  an  order-of-magnitude  comparison  of  these  terms,  the  LPEs  may  be  used 
directly.  Also,  since  Rm  >  R»^  only  Rm  will  be  used  for  this  comparison.  Using  Eqs  (2), 
(4),  and  (19): 

dR^  Zn^J2Rl  ...  /,,, 

- - IrrSmiCOSl  (51) 

9^  2(l-e2)^/2 

^  (52) 

de  2  J 

dQ  3nJ2Rg 

dt  2a2(l-e2)^  ^  ’ 

^  £  sin^  A  +  cos^  i  (54) 

dt  2a2  (1  -  e2)2  V  2  /  ^  2a^  (1  -  ^ 

Considering  only  the  magnitudes  of  the  rates  of  change  of  fl  and  w,  the  following  approx- 


cos^  i 


imations  may  be  obtained: 


dt  ^  dt  ^  aV2(i  _  e2)2 
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Therefore,  averaging  the  effects  of  the  Sun’s  mean  motion  n,  is  valid  when 

n,  > 


a7/2(i  _e2)2 

The  parameters  in  Eq  (56)  have  the  following  approximate  values  (12): 


(56) 


Re 

fJ-m 

n. 


=  0.0019604 
=  3395  km 
=  42828  km^  / sec^ 

=  1.058  X  10"'^  sec~'^ 


The  motivation  for  using  high  eccentricity  orbits  is  to  reduce  orbit  insertion  costs  while 
allowing  for  a  satellite  mission  which  requires  close  approaches.  To  meet  this  requirement, 
consider  the  periapse  radius  Vp  to  be  in  the  range  4000  km  <  Tp  <  7000  km  (it  will  be 
shown  surface  impacts  due  to  the  variation,  of  eccentriciy  do  not  allow  the  examination  of 
periapse  radii  smaller  than  4000km).  Using  the  relation  a  =  rp/(l-e),  averaging  the  solar 
central  angle  requires 


n,  > 


(57) 


Figure  9  presents  a  plot  of  <3  vs  e  for  various  values  of  periapse  radius,  along  with  the 
plot  of  n,  for  comparison.  The  figure  indicates  averaging  the  effects  of  the  Sun  is  valid  for 
high-eccentricity  orbits,  with  the  lower  bound  on  eccentricity  increasing  with  decreasing 
periapse  radius. 
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Eccentricity 


Figure  9. 


Approximate  Magnitude  (Q)  of  u>  and  n  vs  e  for  Various  Values  of  Periapse 
Radius 


Averaging  the  motion  of  the  Sun  within  the  range  determined  above  produces  a 
disturbing  function  which  no  longer  contains  time- varying  solar  elements.  The  Hamiltonian 
for  this  system  may  be  formed  (10,  8,  6),  allowing  for  a  czmonical  transformation  into  a 
new  set  of  variables  in  which  the  new  disturbing  function  is  expandable  in  powers  of  J2- 
Higher  order  terms  are  then  ignored,  resulting  in  a  disturbing  function  containing  critical 
divisors  (7),  all  of  which  are  functions  only  of  inclination.  The  resulting  critical  inclinations 
which  lie  between  0°  and  90®  are  given  in  the  table  below. 
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Table  1.  Critical  Inclinations  from  Analytical  Approach 


ic  (degrees) 


46.4 


56.1 


63.4 


69.0 


73.2 


The  numerical  approach  which  follows  does  not  invoke  the  simplihcations  of  solar 
averaging  and  includes  all  terms  related  to  J2.  The  results  of  the  numerical  integration 
will  be  shown  to  depart  from  the  analytical  results  given  in  the  table  above. 
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III.  Numerical  Development  and  Results 


Numerical  Solution  of  Lagrange ’s  Planetary  Equations 

Using  the  disturbing  function  derived  in  the  last  chapter,  the  LPEs  may  now  be 
numerically  integrated.  Integration  was  accomplished  using  a  multi-step,  variable  step  size 
and  variable  order  scheme  (see  [8]  and  [11]).  To  avoid  numerical  difficulties  which  arise 
when  the  classical  variable  set  (a,e,t,u;,fl,  Af)  is  used,  the  following  transformations  are 
utilized: 

h  =  esinw  (58) 

k  =  ecosu  (59) 

and 

An  =  +  (n  -  <A)  (60) 

•30 

where  <f>  is  the  Mars  equivalent  of  the  Greenwich  hour  angle  and  So  approximates  the 
ratio  of  the  number  of  orbits  per  planet  revolution.  The  variables  h  and  k  axe  introduced 
in  order  to  avoid  the  singularity  at  zero  eccentricity.  The  variable  An  ,  which  is  called 
the  stroboscopic  mean  node  [6],  represents  the  mean  satellite  position.  Introduction  of  An 
allows  for  averaging  the  mean  anomaly  while  retaining  the  effects  of  resonances  with  tesseral 
harmonics.  If  the  position  of  the  satellite  is  commensurate  with  one  or  more  tesserals,  these 
tesserals  will  have  long-term  effects  on  the  orbital  elements  (commensurability  occurs  when 
So  is  a  ratio  of  two  integers). 

Prior  to  casting  the  LPEls  in  terms  of  the  new  variable  set  (o,  h,  t,  fc,  H,  An),  note  the 
disturbing  function  of  interest  does  not  include  the  mean  anomaly,  which  was  averaged 
out.  Therefore,  the  semi-major  axis  a  remains  constant  {da/dt  =  0).  Also,  since  J2  is 
the  only  harmonic  considered,  the  stroboscopic  mean  node  is  not  needed.  The  six  original 
LPEs  (Eqs.  2  to  7)  reduce  to  four  (Elqs.  2,  3,  4,  and  6).  Transforming  these  four  into  the 
new  variable  set  requires  the  following  derivatives: 


de 


% 


dh  de  dk  de 
h dR  k dR 

e  dh  ^  e  dk 


(62) 


dh  dh  de  dh  dw  hde  du> 

dt  de  dt  ^  duj  dt  e  dt"^  dt 


b  dR  k  cos  i  dR 
na^  dk  nab  sin  i  di 


dk  dk  de  dk  dw  k  de  ^dw 

dt  de  dt  ^  du  dt  e  dt  dt 


b  dR  cos  i  dR 
na^  dh  nab  sin  i  di 


cost  dR  b  /'hdR  kdRY 
nab  sin  i  di  na^e  \e  dh  e  dk  J  _ 

(63) 


cost  dR  b  /hdR  kdRY 
nab  sin  i  di  na^e  \e  dh  e  dk  J  _ 

(64) 


In  terms  of  the  new  variable  set,  the  pertinent  LPEs  become 


dil 

1  dR 

dt 

nab  sin  i  di 

di 

cot  i  /  dR 

1  dR 

dt 

nab  V  dh 

^dk) 

nab  sin  i  d^l 

dh 

b  dR  k  cot  i  dR 

dt 

na^  dk 

nab  di 

dk 

b  dR 

h  cot  i  dR 

dt 

na^  dh 

nab  di 

where  6  =  a(l  —  is  the  semi-minor  axis  as  before,  and 

e  =  /lih^  +  k^) 


(65) 

(66) 

(67) 

(68) 


(69) 


The  complete  disturbing  function  in  terms  of  the  new  variables  is 


R  = 


+ 

-I- 

+ 


_ _ 

2a3(l  _/i2  _  fc2)3/2 


(' 


^  {y  [ciC^_A  -  (CiCr,5n-A  -  SlS^)^]  {k^  ~  h^) 

30  ^CiSiSiCn-A  -  C'tC.Cn-A'S'n-A]  hk 

(2  -f  3h^  +  3k^)  [clCl_^  +  {CiCiSa-A  -  5l5.)2]  -  l]  | 


(70) 
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The  required  partial  derivatives  are 


dR 

dk 


dR 

dh 


2a3(l-/i2-A:2)5/2 


+  {l5  [ciCS-A  -  (CiCAS„.A  -  5Lii)’l  * 

+  30  [c^SiSiCo-A  -  ClCiCii-KSo-K^  h 
+  6k  [c^cLa  +  (CiCLSa-A  -  5i5.)"]  -  l]  } 


Z^imhRlh 


2a^(l-h^  -k^f^ 


(l-2sln^) 


+  ^{-i5[c£C^_^-(C.Cz,5n_A-5L50' 

+  30  [C£,5L-^i<^n-A  -  CiCiCu-KS^.t^  k 

+  6h  [clCl_^  +  {CiCLS(i-K  -  Sl5.)"]  -  l]  I 


(71) 


(72) 


di 


3/^  Jj-^e 


4a3  (1  _  /r2  _  fc2)3/2 


sin  2i 


ti»a 

4rf 


{l5(CiCLSn_A  -  SiSi)(SiCLSn.A  +  SlQ)  (k^  -  h^) 
+  30  [Cl5i,C, Cn-A  +  C’i'^iC'n-A-^n-A] 

-  3  (2  +  3/1^  +  3fc2)  {CiCcSu-A  -  SiSi)  {SiCiSn-A  +  ^lQ)} 


(73) 


^  I — Y  [C'i52(n_A)  +  2(CiCL5n_A  -  ■5x,5,)C,Cf,Cn_A]  (fc^  -  h^) 

—  30  ^C'£,5£,SjS’n-A  +  C'2CiC2(n_A)j  hk 

-  ^  (2  +  3/1^  +  3k^)  [C£52(n-A)  -  2  (CiCiSu-A  -  5l5.)  C.CiCn_A] }  (74) 

The  details  of  employing  the  numerical  integrator  to  locate  critical  inclinations  is 
described  in  the  next  section. 
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Locating  Critical  Inclinations 

Figures  1  through  4  illustrated  the  phenomenon  of  large  amplitude  variations  in  ec¬ 
centricity  and  inclination.  These  large  variations  occur  under  specific  conditions,  namely 
certain  critical  inclinations  for  high  eccentricity  orbits.  The  last  section  of  Chapter  II  sum¬ 
marized  an  analytical  approach  for  determining  the  location  of  these  critical  inclinations, 
which  were  found  to  be  independent  of  eccentricity  (i.e.,  constant).  The  resultant  criti¬ 
cal  inclinations  also  had  the  property  of  producing  a  local  maximum  in  the  variations  of 
both  eccentricity  and  inclination.  The  intent  of  the  numerical  approach  is  to  accomplish 
a  discrete  numerical  search  for  the  critical  inclinations  throughout  a  range  of  interest.  For 
the  purpose  of  this  search,  a  critical  inclination  is  defined  as  the  inclination  at  orbit  inser¬ 
tion  which  results  in  the  largest  local  variation  in  eccentricity  or  inclination  (i.e,  a  local 
maximum).  No  assumption  is  made  concerning  correlations  between  the  behavior  of  the 
eccentricity  and  the  behavior  of  inclination,  therefore  a  separate  search  is  made  for  each. 

Search  Method. 

The  numerical  integrator  produces  an  output  data  file  consisting  of  the  values  of  the 
orbited  elements  at  time  increments  determined  by  the  user.  The  behavior  of  the  elements 
may  be  observed  by  plotting  the  data.  However,  the  magnitude  of  the  search  demands  the 
use  of  a  numerical  technique  for  locating  the  local  maximums  in  the  variations,  thereby 
identifying  the  critical  inclinations.  The  method  employed  is  to  fit  a  linear  function,  in 
the  least  squares  sense  [9:23-47],  to  the  graph  of  eccentricity  versus  time  and  the  graph  of 
inclination  verses  time.  The  standard  deviation  of  the  residuals  then  becomes  the  tool  for 
locating  critical  inclinations.  The  larger  the  variation  in  the  behavior  of  these  elements, 
the  greater  the  standard  deviation  will  be  for  the  linear  fit.  By  looping  through  the 
range  of  orbit  insertion  inclinations  while  keeping  other  initial  conditions  the  same,  critical 
inclinations  for  a  given  eccentricity  are  those  which  produce  a  local  maximum  in  the  plot  of 
standard  deviation  verses  inclination.  An  outer  loop  which  increments  the  orbit  insertion 
eccentricity  completes  the  search,  allowing  for  a  three-  dimensional  surface  to  be  plotted, 
where  eccentricity  and  incUnation  are  the  independent  variables,  and  standard  deviation 
is  the  dependent  variable.  The  ‘peaks’  in  this  surface  are  the  points  of  interest. 
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Search  Parameter  Values  and  Initial  Conditions. 

The  search  method  described  above  would  be  best  accomplished  with  very  small 
increments  in  eccentricity  and  inclination,  combined  with  a  sufficiently  long  propagation 
time.  The  plot  in  Figure  3,  showing  data  which  spans  a  period  of  seventy  Earth  years 
about  Mars,  does  not  yet  exhibit  a  repeatable  pattern.  However,  computer  time  and 
storage  constraints  did  not  allow  for  large  data  collection,  therefore  a  smaller  sample  was 
used.  This  is  not  a  serious  constraint  since  large  variations  near  critical  inclinations  become 
evident  very  soon  after  orbit  insertion.  Also,  an  orbit  period  of  five  or  ten  years  may  be 
of  more  immediate  concern. 

The  table  below  lists  the  initial  conditions  and  search  parameter  increments  used  to 
perform  the  numerical  search.  Choosing  the  peri  apse  radius  Tp  determines  the  value  of  the 
semi-major  axis  by  a  =  rp/(l  —  e).  For  a  complete  listing  of  the  program  input  file,  see 
appendix  B. 

Table  2.  Search  Increment  Values  and  Initial  Conditions 


w  at  insertion 

0 

n  at  insertion 

0 

M  at  insertion 

0 

e  at  insertion 

incremented  from  0.40  to  0.90  by  0.02 

i  at  insertion 

incremented  from  0.25  to  90.0  by  0.25  (degrees) 

Tp  at  insertion 

incremented  from  4500  km  to  7000  km  by  500  km 

orbit  lifetime 

10  Earth  years 

epoch  at  insertion 

7  January  2001 

data  output  increment 

100  days 

The  effect  of  using  a  different  epoch  at  orbit  insertion  was  tested  and  found  to  produce  no 
noticeable  change  in  the  results  for  the  search  range  considered. 


Critical  Inclinations  for  Eccentricity 


The  inner  loop  of  the  search  procedure  described  in  the  last  section  consisted  of 
incrementing  through  the  range  of  inclinations  for  a  given  eccentricity.  For  each  of  these 
insertion  eccentricities,  a  plot  of  the  standard  deviation  of  the  residuals  in  eccentricity 
(SDE)  verses  inclination  may  be  generated.  These  plots  are  too  numerous  to  include  them 
all  here,  however,  all  such  plots  for  a  periapse  radius  of  6500A:m  can  be  found  in  appendix  D. 


Three-Dimensional  Surface  Plots. 

Combining  all  the  plots  of  SDE  vs  Inclination  for  a  given  periapse  radius  allows  for 
the  generation  of  a  three-dimensional  plot,  where  the  third  axis  is  the  range  of  eccentricities 
(these  plots  are  presented  in  Figures  10  through  16).  These  surfaces  indicate  the  locations 
of  the  critical  inclinations,  which  vary  with  eccentricity. 


Critical  inclinations  at  eccentricities  near  0.4  and  inclinations  near  0.25  degrees  vary 
rapidly  with  eccentricity,  generating  a  rough  surface  topography.  This  rough  area  moves 
out  of  the  search  range  as  periapse  radius  increases  and  is  totally  absent  at  periapse  radii 
greater  than  6000km.  Subsequent  analysis  will  therefore  focus  on  the  remaining  smoother 
areas  where  pronounced  topography  is  present  throughout  the  range  of  interest. 


Some  of  the  surfaces  do  not  contain  the  complete  range  of  data  (i.e.,  eccentricities 
between  0.4  and  0.9).  Where  this  data  is  not  present,  large  variations  in  the  eccentricity 
caused  the  orbit  to  impact  before  the  full  set  of  data  was  collected.  These  computer  runs 
were  rejected  (see  appendix  C  for  further  discussion). 
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Two-Dimensional  Plots. 


In  order  to  reduce  each  three-dimensional  surface  to  a  two-dimensional  plot  of  critical 
inclination  verses  eccentricity,  the  following  procedure  was  employed.  For  each  eccentricity 
considered,  local  maximums  were  identified  by  observing  the  slope  of  the  line  connecting 
adjacent  data  points  in  the  plot  of  SDE  versus  inclination.  A  data  point  was  considered 
a  maximum  (critical)  when  the  slope  of  the  line  transitioned  from  positive  to  negative. 
Figures  10  through  16  indicate  the  presence  of  five  significant  critical  inclination  curves, 'as 
well  as  other  smaller  local  maximums.  These  smaller  maximums  were  filtered  out,  retaining 
maximums  with  the  five  greatest  values  of  SDE.  The  resultant  data  is  given  in  tabular  form 
in  appendix  E,  and  the  plots  of  critical  inclinations  verses  eccentricity  are  given  in  Figures 
17  through  22.  Insufficient  data  was  collected  at  periapse  radius  of  4000  km  to  produce  a 
meaningful  two-dimensional  plot,  therefore  only  periapse  radii  of  4500  km  and  above  are 
presented. 
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Figure  18.  Critical  Inclination  vs  Eccentricity;  Periapse  Iladius=5000  km 
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Figure  19.  Critical  Inclination  vs  Eccentricity:  Periapse  RaKlius=5500  km 
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Critical  Inclinations  for  Inclination 

The  techniques  for  locating  critical  inclinations  in  eccentricity  were  also  applied  to 
inclination.  Incrementing  through  the  range  of  inclinations  for  a  given  eccentricity  produces 
a  plot  of  the  standard  deviation  of  the  residuals  in  inclination  (SDI)  verses  inclination  (see 
appendix  D).  As  before,  combining  the  data  for  the  full  range  of  eccentricities  allows  for 
the  creation  of  a  three-dimensional  surface  plot  for  each  periapse  radius  considered.  These 
surface  plots  are  given  in  Figures  23  through  29. 

Unlike  the  surfaces  for  eccentricity,  the  surfaces  for  inclination  indicate  fewer  distinct 
critical  locations.  The  presence  of  one  or  two  well  pronounced  critical  locations  in  the 
surface  topography  at  high  eccentricities  forces  other  possible  critical  locations  to  appear 
flat.  Attempting  to  reduce  the  surfaces  for  inclination  to  two  dimensions  as  was  done  for 
eccentricity  did  not  produce  meaningful  results. 

Although  the  attempt  to  produce  a  two-dimensional  view  of  the  critical  values  for 
inclination  was  not  successful,  overlapping  the  plots  of  SDE  and  SDI  verses  inclination 
for  each  eccentricity  allows  for  a  comparison  of  the  location  of  the  local  maximums.  Ap¬ 
pendix  D  presents  these  plots  for  a  periapse  radius  of  6500km.  There  is  no  indication 
of  any  correlation  between  the  location  of  critical  inclinations  for  eccentricity  and  critical 
inclinations  for  inclination.  However,  recall  these  results  apply  to  a  ten  Earth-year  orbit 
lifetime,  which  does  not  capture  a  full  period  of  the  large  amplitude  variations.  A  longer 
orbit  lifetime  may  produce  different  results. 
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Modeling  the  Locations  of  Critical  Inclinations  in  Eccentricity 

The  locations  of  critical  inclinations  for  each  periapse  radius  considered  were  pre¬ 
sented  in  Figures  17  through  22.  The  shapes  of  the  curves  are  similar  between  periapse 
radii,  but  shifted  slightly  from  one  periapse  radius  to  another.  This  s'  "*;ion  will  take  a 
closer  look  at  the  dependence  of  these  critical  inclination  curves  on  eccentricity  and  peri¬ 
apse  radius,  and  a  curve  fit  will  be  applied  to  three  of  the  curves. 

The  Variation  of  Critical  Inclination  Curves  with  Periapse  Radius. 

To  better  view  the  dependence  of  the  critical  inclination  curves  on  periapse  radius, 
Figures  30  through  35  present  plots  of  the  overlap  of  all  similar  curves  (i.e.,  all  curves 
labeled  ‘a’  appear  in  Figure  30,  all  curves  labeled  ‘b’  in  Figure  31,  and  so  forth).  Each  set 
of  curves  is  similar  in  shape  but  is  shifted  toward  lower  eccentricities  as  periapse  radius 
increases.  Two  other  features  are  also  present. 

The  first  is  the  tendency  of  each  group  of  curves  to  converge  at  lower  eccentricities. 
By  inspection,  curves  a,  b,  and  f  converge  to  approximately  67.5,  64.5,  and  61  degrees 
respectively,  and  curves  e  appear  to  converge  to  approximately  54  degrees.  Curves  ‘c’  and 
‘d’  do  not  extend  far  enough  into  the  lower  eccentricities  to  estimate  their  convergence 
points.  Comparing  these  convergence  inclinations  directly  with  the  critical  inclinations 
determined  by  analytical  means  may  not  be  justified  since  the  numerical  search  did  not 
capture  the  full  period  of  the  variations.  Nevertheless,  curves  ‘a’  may  be  compared  with 
the  69.0  degree  analytical  result,  curves  ‘b’  with  the  63.4  degree  analytical  result,  and 
curves  ‘e’  with  the  56.1  degree  analytical  result. 

The  second  feature  of  the  critical  inclination  curves  is  the  direction  in  which  the 
curves  diverge  from  the  convergence  points  discussed  above.  Curves  a,  b,  and  c  diverge 
toward  increasing  inclinations;  curves  d,  e,  and  f  diverge  toward  decreasing  inclinations. 
The  dividing  inclination  seems  to  be  near  the  well  known  critical  inclination  produced  by 
J2  alone,  namely  63.4  degrees,  at  which  the  argument  of  pericenter  remains  constant.  The 
direction  of  divergence  may  be  related  to  the  direction  of  the  change  in  the  argument  of 
pericenter. 
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Figure  30.  Critical  Inclination  vs  Eccentricity:  Curve  a 
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Figure  31.  Critical  Inclination  vs  Eccentricity:  Curve  b 
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Figure  32.  Critical  Inclination  vs  Eccentricity;  Curve  c 
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Figure  33.  Critical  Inclination  vs  Eccentricity:  Curve  d 
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Figure  34.  Critical  Inclination  vs  Eccentricity:  Curve  e 
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Figure  35.  Critical  Inclination  vs  Eccentricity:  Curve  f 


Curve  Fitting. 

Of  the  six  critical  inclination  curves  plotted,  curves  a,  b,  and  e  best  demonstrate  the 
convergence  behavior  discussed  above.  For  curves  a  and  b,  as  eccentricity  increases  from 
0,4,  the  critical  inclination  remains  constant  until  a  specific  eccentricity  is  reached,  hence¬ 
forth  called  the  departure  eccentricity.  Beyond  the  departure  eccentricities,  the  critical 
inclination  curves  are  no  longer  constant.  Curve  e  exhibits  the  same  behavior,  however  the 
departure  eccentricities  appear  to  occur  below  the  range  of  data  examined. 

In  order  to  apply  multiple  regression  to  curves  a,  b,  and  e,  each  curve  must  be 
divided  into  two  sections.  The  section  which  lies  below  the  departure  eccentricity  is  simply 
a  straight  line  at  constant  inclination.  The  section  above  the  departure  eccentricity  is  the 
section  to  be  fitted  (for  curve  e,  the  entire  range  of  data  will  be  assumed  to  lie  above  the 
departure  eccentricity).  The  location  of  the  departure  eccentricities  must  be  determined 
as  a  function  of  periapse  radius.  Since  the  data  for  each  curve  is  discrete,  the  departure 
eccentricities  for  the  data  collected  may  be  approximated  (recall,  the  search  increment  in 
inclination  was  0.25  degrees).  Table  3  below  lists  the  departure  eccentricities  for  curves  a 
and  b  for  each  periapse  radius. 


Table  3.  Departure  Eccentricities  for  Curves  a,  b,  and  e 


Curve  a 

Curve  b 

Periapse  Radius  (km) 

Departure  Eccentricity 

4500 

0.74 

0.78 

5000 

0.68 

0.74 

5500 

0.62 

0.70 

6000 

0.56 

0.64 

6500 

0.50 

0.58 

7000 

0.44 

0.54 

59 


A  linear  fit  to  the  data  above  yields  correlation  coefficients  of  -1.0000  and  -0.997  for  curves 
a  and  b  respectively.  The  resulting  functions  are 


DEa  =  1.280  -  1.20  X  10'‘‘r, 
DEb  =  1.235  -  9.94  x  lO'^r, 


(75) 

(76) 


where  DEa  and  DEb  are  the  respective  departure  eccentricities  for  curves  a  and  b,  and 
Tp  is  in  kilometers. 

Fitting  the  critical  inclination  curves  above  the  departure  eccentricities  to  a  curve 
linear  in  periapse  radius  and  quadratic  in  eccentricity  yields  the  following  functions: 


^ca  “  \ 


^cb  —  \ 


67.5  c  <  DEa 
-33.9  +  0.01204rp  -f  89.4036^  e  >  DEa 

64.5  e  <  DEb 
0.0068443rp  -|-  56.41  le^  e  >  DEb 

»ce  =  I  107.72  -  0.0067837rp  -  49.59e2  e  >  0.40 

where  ica  ,  icb  >  and  ice  are  the  respective  critical  inclination  functions  for  curves  a,  b, 
and  e.  The  coefficients  of  multiple  determination  (R^)  for  these  fits  axe  0.95,  0.94,  and 
0.86  respectively.  The  curve  fits  were  accomplished  using  the  GLM  (general  linear  model) 
procedure  of  the  SAS  (Statistical  Analysis  System)  computer  package  [10]. 
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IV.  Conclusions  and  Recommendations 
Summary  and  Conclusions 

The  purpose  of  this  thesis  was  to  locate  critical  inclinations  in  long  term  high  eccen¬ 
tricity  orbits  about  Mars  using  numerical  methods.  The  approach  consisted  of  applying  a 
linear  lea.st  squares  fit  to  the  graphs  of  eccentricity  versus  time  and  inclination  versus  time 
for  a  range  of  eccentricities  and  inclinations.  The  residuals  from  the  linear  fits  served  as 
the  search  parameter  for  identifying  the  critical  inclinations,  which  appeared  in  the  form 
of  local  maximums  in  the  three-dimensional  surface  plots  of  the  standard  deviation  in  the 
residuals  versus  e  and  i.  The  following  results  were  obtained: 

1.  The  numerical  approach  clearly  identified  the  locations  of  some  of  the  critical  in¬ 
clinations  in  eccentricity  within  the  range  searched.  The  presence  of  other  critical 
inclinations  were  indicated,  however  the  range  and  resolution  of  the  search  were  not 
extensive  enough  for  close  examination  of  these  other  locations. 

2.  The  locations  of  critical  values  for  eccentricity  were  found  to  vary  with  eccentricity 
and  periapse  radius.  Six  distinct  critical  location  curves  were  identified,  three  of 
which  were  curve  fit  to  functions  linear  in  periapse  radius  and  quadratic  in  eccentric¬ 
ity. 

3.  These  three  curves  converge  to  three  critical  inclinations  determined  by  analytical 
results  to  within  the  search  increments  used.  The  remaining  curves  did  not  include 
enough  data  for  direct  comparison  with  analytical  results. 

4.  Critical  values  for  inclination  were  not  clearly  identified.  Three  dimensional  surface 
plots  of  the  standard  deviation  in  the  residuals  for  inclination  produced  one  or  two 
dominant  local  maximums. 

5.  For  the  range  and  resolution  of  data  collected,  no  evidence  of  correlation  between 
critical  values  for  eccentricity  and  critical  values  for  inclination  was  found. 

Based  on  these  results,  a  high  eccentricity  orbit  about  Mars  may  be  utilized  to 
meet  specific  satellite  mission  requirements,  or  to  reduce  orbit  insertion  costs,  provided 
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care  is  taken  to  avoid  critical  inclinations  which  produce  unacceptable  variations.  This 
investigation  demonstrated  a  numerical  procedure  for  identifying  the  conditions  for  which 
such  variations  occur. 

Recommendations  for  Further  Study 

The  presence  of  constraints  within  any  form  of  research  automatically  provides  op¬ 
portunities  for  further  study  when  the  contraints  are  removed.  The  work  presented  here 
could  be  repeated  using  much  finer  search  increments,  and  more  importantly,  longer  orbit 
lifetimes.  Capturing  at  least  one  full  period  of  the  variations  in  eccentricity  and  inclination 
may  provide  for  a  much  more  detailed  analysis,  as  well  as  a  more  complete  comparison 
with  analytical  results.  Under  these  more  favorable  conditions,  the  following  additional 
studies  might  also  be  accomplished: 

1.  Attempt  to  determine  low  eccentricity  convergence  points  for  all  critical  inclination 
curves  identified. 

2.  Determine  not  only  the  location  of  the  critical  inclinations  but  also  the  magnitudes 
and  periods  of  the  variations. 

3.  Extend  the  search  range  to  180  degrees  inclination.  The  work  accomplished  in  this 
investigation  produced  no  evidence  of  symmetry  about  90  degrees  inclination. 

4.  Include  the  effects  of  additional  harmonics  in  the  Mars  gravity  potential. 

5.  Investigate  the  effects  of  varying  the  arguments  of  periapse  and  ascending  node  at 
orbit  insertion. 
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Appendix  A.  Hansen’s  Coefficients 


For  a  complete  derivation  of  Hansen’s  Coefficients,  see  [2].  This  appendix  provides 
a  brief  definition,  as  well  as  the  derivation  of  Eqs  (14),  (17),  and  (18),  and  Eqs  (45),  (46), 
and  (47). 

Definition  of  Hansen’s  Coefficient 
Let 

X  =  exp  (jf) ,  y  =  exp  (jE) ,  2  =  exp  (jAf ) 

where  /  is  the  true  anomaly,  E  is  the  eccentric  anomaly,  M  is  the  mean  anomoly,  and 
j  =  y/—l-  Hansen’s  Coefficients  are  then  defined 

For  p  =  0; 

xr  =  (-!)"■'  f  "  ^  )  C  ,  I2L1I;  ;  |m|  +  1  ,  e=)  (78) 

where  F  is  the  hypergeometric  function  [1:272-277]  defined  as 

The  abbreviation  (a)„,  which  is  called  the  Pochhammer  symbol,  has  the  following  proper¬ 


ties: 


(c)o  =  1 

(l)n  =  n! 

(a)„  =  a(a-|-  l)(a -f- 2)  •  • -(a n  -  1)  n  =  1,2,3,... 
The  Pochhammer  symbol  and  the  binomial  coefficient  are  related  by 


(a)„  =  (-!)’*«! 


—a 


n 


The  following  property  of  Hansen’s  Coefficients  will  be  used  below; 

yTHjTl  vTU,— tl 

Aq  —  Aq 


(79) 
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Derivations 


Derivation  of  Eq  (14); 


The  last  expression  is  simply  the  binomial  expansion  for  (1  -  so 


Derivation  of  Eq  (17); 

Using  the  identity 

sin  2f=-^[  exp(j2/)  -  exp(->2/)  ] 

then 

=  hi^riz)  '“^->2/)" 


because  of  Eq  (79). 
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Derivation  of  Eq  (46) 


27r 


2r 


2 

sin  2  /  dM 


=  i  [irr"  -  i;C  {^y^M-nndM 

=  i  -  xi--'‘)  =  0  (88) 

Derivation  of  Eq  (47) 
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Appendix  B.  Sample  Program  Input  File 


This  appendix  includes  the  program  input  file  for  the  numerical  integrator  (this 
information  is  included  to  provide  consistency  should  further  analysis  be  performed).  Only 
A  (a),  E  (e),  and  I  (i)  were  changed  for  each  computer  run.  All  other  parameters  remained 
constant.  For  an  explanation  of  the  parameters  listed,  see  [8]. 


70000. DO 


L 

M 

ISUN 

IMOON 

lEPHEM 

IDRAG 

ISRP 

I  ORB 

IPRINT 

IPLOT 

NP 

NQ 

IQMAX 

NK 

ICASE 

LSUN 

NSUN(l) 

(2) 

(3) 

(4) 

(5) 

(6) 

(7) 

(8) 

(9) 

LMOON 

NMOON(l) 

(2) 

(3) 

(4) 

(5) 

(6) 
(7) 
(6) 
(9) 

NSEGl 
NSEG2 
ORB(l),  A 
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0.9D0 

(2). 

40.75D0 

(3). 

O.DO 

(4), 

O.DO 

(5), 

O.DO 

(6). 

l.D-6 

RELERR 

l.D-6 

ABSERR 

100. DO 

STEP 

19911007. DO 

TINT(l) 

O.DO 

(2) 

200 11007. DO 

TFIN(l) 

O.DO 

(2) 

19911007. DO 

TREF(l) 

O.DO 

(2) 

4.2828287D4 

GE 

3397. 2D0 

RE 

4.061249803D-3 

RATE 

333. 5597 IDO 

PM 

.1017D0 

ELLIP 

3487. 2D0 

RATM 

6.D-4 

RDENS 

361. DO 

RHT 

36.  DO 

SHT 

1000. DO 

ALTMAX 

lO.D-6 

AREAD 

20.D-6 

AREAS 

1000. DO 

SCHASS 

2. DO 

CDRAG 

1.95D-3 

CSRP 

.13271244D12 

GS 

227.9410D6 

ES(1) 

93.39697D-3 

ES(2) 

25.191153D0 

ES(3) 

O.DO 

ES(4) 

-109.0S06D0 

ES(5) 

171.60476D0 

ES(6) 

6.065196184D-6 

ES(7) 

O.DO 

GH 

O.DO 

EM(1) 

O.DO 

EM(2) 

O.DO 

EM(3) 

O.DO 

EM(4) 

O.DO 

EM(5) 

O.DO 

EM(6) 

O.DO 

EM(7) 
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-0 . 1960387250000000E-02 


0 . OOOOOOOOOOOOOOOOE+00 
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Appendix  C.  Linear  Least  Squares  Fit  for  Eccentricity  and  Inclination 


Linear  fits  to  the  curves  of  eccentricity  versus  time  and  inclination  versus  time  allowed 
for  using  the  standard  deviation  of  the  residuals  as  the  search  parameter  for  locating  critical 
inclinations.  In  the  standard  notation,  x  is  the  independent  variable,  y  is  the  dependent 
variable,  and  s  is  the  residual  (13:23-47),  so  that 


Vi  =  /?o  +  0ix,  +  Ci 


(90) 


where  the  subscript  i  refers  to  the  ith  observation  or  data  point.  The  coefRcients  are 
determined  bv 


Pi  —  s2 

X  n 


(91) 

(92) 


where  n  is  the  number  of  data  points.  The  error  mean  square  (MSE)  of  the  residuals  is 
given  bv 

i:iyi-0o-0iXif  _ 


MSE  = 


(93) 


n  —  2  n  —  2 

The  estimator  of  the  standard  deviation  of  the  residuals  is  simply  the  positive  square  root 
of  MSE. 


A  FORTRAN  program  was  written  to  calculate  the  standard  deviations.  The  pro¬ 
gram  reads  data  from  the  numerical  integrator  output  file,  calculates  the  standard  deviation 
of  the  residuals  to  the  linear  fit  in  eccentricity  and  inclination  for  each  orbit  run,  then  out¬ 
puts  the  standard  deviations  to  a  new  file.  Three-dimensional  surface  plots  were  created 
using  the  data  from  this  new  file. 

The  program  includes  a  provision  for  computer  runs  which  ended  prematurely.  If  the 
radius  of  periapse  became  equal  to  or  less  than  the  planet  radius  due  to  a  large  variation  in 
eccentricity,  an  impact  occured  and  the  run  was  terminated.  In  order  to  use  the  standard 
deviation  as  a  search  parameter  for  critical  inclinations,  each  run  must  have  the  same  (or 
nearly  the  same)  number  of  data  points.  Significantly  shorter  runs  were  rejected,  runs 
which  were  only  a  few  days  shorter  were  kept.  The  shorter  runs  appear  as  discontinuities 
at  the  edges  of  the  surface  plots. 
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Appendix  D.  Plots  of  SDE  and  SDI  versus  Inclination  for  rp  =  6500km 


This  appendix  includes  all  the  plots  of  Standard  Deviation  in  Eccentricity  (SDE) 
and  Standard  Deviation  in  Inclination  (SDI)  versus  Inclination  for  a  periapse  radius  of 
6500km.  The  plots  of  SDE  and  SDI  are  overlayed  for  the  purpose  of  comparison.  Each 
vertical  axis  is  relatively  scaled. 
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Figure  42.  SDE  and  SDI  (dashed)  vs  Inclination 
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Figure  43.  SDE  and  SDI  (dashed)  vs  Inclination 
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Figure  45.  SDE  and  SDI  (dashed)  vs  Inclination 
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Figure  46.  SDE  and  SDI  (dashed)  vs  Inclination 
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Figure  52.  SDE  and  SDI  (dashed)  vs  Inclination 
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Figure  53.  SDE  and  SDI  (dashed)  vs  Inclination 


Eccentric! t  y-0 . 84 


SDE 
0 . 0052 
0 . 0050 
0048 
0046 
0044 
,0042  -1 
,  0040 
0.0038  - 
0.0036  - 
.0034  - 
.  0032  - 
,  0030 
.0028  - 
.0026  - 
.0024  - 
0.0022  - 
0.0020  ■ 
0.0018  - 
0.0016  - 


SOI 


0. 
0  . 
0  , 
0. 
0, 


0. 
0. 
0. 
0  . 
0  . 
0. 


30  40 

Inclination  (degrees) 

Figure  58.  SDE  and  SDI  (dashed)  vs  Inclination 
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Figure  59.  SDE  and  SDI  (dashed)  vs  Inclination 
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Figure  61.  SDE  and  SDI  (dashed)  vs  Inclination 
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Appendix  E. 


Tables  for  Critical  Inclinations  i: 


Eccentricity 
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Table  4.  Critical  Inclinations  for  Eccentricity:  Curve  a 


Radius  of  Periapse  (km) 


4500  5000  5500  6000  6500  7000 


Critical  Inclinations  (degrees) 


67.75  68.00  68.25 


67.50  68.00  68.00  68.50 


67.50  67.75  68.25  68.75 


67.75  67.75  68.25  69.25 


67.75  68.00  68.50  70.00 


67.50  67.75  68.25  68.75  71.75 


67.50  67.75  68.25  69.25  74.00 


67.75  68.00  68.50  70.25  75.75 


68.00  68.25  68.75  72.00  77.50 


67.50  67.75  68.25  69.25  74.25  79.50 


67.50  67.75  68.50  70.25  76.25  81.75 


67.75  68.00  68.75  72.50  78.25  85.25 


68.00  68.25  69.25  75.00  80.50 


67.75  68.50  70.50  77.00  83.75 


68.00  68.75  73.25  79.50 


68.25  69.50  75.75  82.50 


68.50  71.25  78.25 


68.75  74.50  81.50 


69.75  77.25  87.50 


73.00  80.75 


76.50  87.75 


80.50 


Table  5.  Critical  Inclinations  for  Eccentricity:  Curve  b 


Radius  of  Periapse  (km) 


4500  5000  5500  6000  6500  7000 


Critical  Inclinations  (degrees) 


65.00  64.75 


65.00  64.75 


64.75  64.50 


64.75  64.50 


65.00  64.75  64.75 


65.00  64.50  64.75 


64.75  64.75  64.50  65.00 


65.00  64.75  64.75  65.75 


65.00  64.75  64.75  66.50 


64.75  64.75  65.25  67.50 


64.75  64.75  65.75  68.75 
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65.00  64.75  65.00  66.75  69.75 


65.00  64.75  65.50  68.00  70,75 


64.50  64.75  64.75  66.25  69.50  71.75 


65.00  64.75  65.00  67.50  70.50  73.00 


65.00  64,75  65.75  69.00  71.75  74.00 


64.75  65.00  66.75  70.25  72.75  75.50 


64.75  65.25  68.50  71.50  74.00  77.25 


64.75  66.25  70.00  72.75  75.75  79.50 


65.25  68.00  71.50  74.25  77.75  84.25 


66.00  70.00  73.00  76.25  81.50 


67.75  71.50  74.75  79.25 
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Table  6.  Critical  Inclinations  for  Eccentricity:  Curve  c 


Radius  of  Periapse  (km) 

Critical  Inclinations  (degrees) 

e 

4500 

5000 

5500 

6000 

6500 

7000 

0.52 

- 

- 

- 

- 

- 

60.25 

0.54 

- 

- 

- 

- 

- 

60.25 

0.56 

- 

- 

' 

- 

- 

60.50 

0.58 

- 

- 

- 

- 

- 

60.50 

0.60 

- 

- 

- 

- 

- 

60.75 

0.62 

- 

- 

- 

- 

60.50 

61.25 

0.64 

- 

- 

- 

- 

60.75 

61.50 

0.66 

- 

- 

- 

60.50 

61.25 

62.00 

0.68 

■ 

- 

- 

60.75 

61.50 

62.25 

0.70 

- 

- 

- 

61.00 

62.00 

62.50 

0.72 

- 

- 

60.75 

61.50 

62.25 

62.50 

m 

- 

- 

61.25 

62.00 

62.50 

62.50 

0.76 

- 

- 

61.75 

62.50 

62.75 

62.75 

0.781 

- 

61.25 

62.25 

63.00 

63.00 

62.75 

m 

- 

62.00 

63.00 

63.25 

63.00 

62.75 

0.82 

62.00 

62.75 

63.50 

63.50 

63.25 

62.75 

0.84 

63.00 

64.00 

63.75 

63.75 

63.25 

62.75 

0.86 

64.00 

64.75 

64.25 

64.00 

63.50 

63.00 

0.88 

65.00 

65.00 

65.00 

64.50 

64.00 

63.50 

0.90 

- 

- 

- 

65.25 

65.00 

64.75 

Table  7.  Critical  Inclinations  for  Eccentricity:  Curve  d 


Table  8.  Critical  Inclinations  for  Eccentricity;  Curve  e 


Radius  of  Periapse  (km) 


4500  5000  5500  6000  6500  7000 


Critical  Inclinations  (degrees) 


54.75 


54.75 


54.50 


53.75  52.25  52.00 


53.50  52.25  51.50 


53.00  52.25  51.25 


54.75 


54.50 


54.50  54.00  52.50  52.00  50.75 


54.25  53.75  52.25  51.75  50.25 


53.25  52.25  51.25  49.50 


52.75  52.00  50.75  49.25 


54.50  54.00  52.25  51.75  50.00  49.25 


54.25  53.75  52.25  51.25  49.50  48.75 


53.25  52.25  50.75  49.25  47.75 


52.50  51.75  50.00  49.25  46.75 


53.75  52.25  51.25  49.50  48.50  45.50 


53.50  52.25  50.75  49.50  47.25  44.25 


52.75  52.00  50.00  49.00  46.00  42.75 


52.25  51.25  49.50  48.00  44.75  41.25 


52.25  50.75  49.75  46.75  43.25  39.25 


52.00  50.00  48.50  45.25  41.50  37.25 


51.25  49.75  47.25  43.50  39.50  34.75 


50.50  49.25  45.75  41.75  37.25  31.75 


49.75  47.75  44.00  39.50  34.25  28.25 


49.75  46.00  41.75  36.75  30.75  24.50 


48.25  44.00  39.25  33.50  26.75  20.50 


46.25  41.50  36.00  29.25  22.25  16.75 


43.75  38.25  31.75  24.50  18.00  14.25 


40.75  34.25  26.50  19.50  15 


m 


Table  9.  Critical  Inclinations  for  Eccentricity:  Curve  f 


Radius  of  Periapse  (km) 


4500  5000  5500  6000  6500  7000 


Critical  Inclinations  (degrees) 


60.25 

60.25 

60.25 

60.25 

60.25 

60.25 

60.25 

60.50 

60.25 

60.00 

60.25 

60.25 

60.25 

60.25 

59.75 

60.25 

60.25 

60.25 

60.25 

59.25 

60.25 

60.25 

60.25 

60.25 

58.75 

60.25 

60.50 

60.25 

1 

60.00 

- 

0.52  60.25  60.25  60.25  60.25  59.50 


0.54  60.25  60.25(60.25  60.25  58.75 


0.56  60.25160.50  1 60.25  60.00 


0.58  60.25  60.50  60.50  59.50 


0.60  60.25  60.50j  60.25  59.00 


0.62  60.50  1 60.50  60.25 


0.64]  60.50  60.50  59.75 


0.66  60.50  60.50  59.00 


0.68  60.50  60.25 


0.70  60.50  60.00 


0.72  60.50  59.25 


0.74  60.50 


0.76  60.00 


Q 


78  59.25 
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